function STA(varargin)
% Spike-Triggered average of membrane potential around the spike
% This returns the modified control worksheet

% Modified 2/10/2000 : changed so that changing rates in the middle of the
% run does not cause subsequent traces to be incorrectly analyzed.
% i.e., all time pointers are recomputed for every trace, depending on the
% current rate...
% Affected also: find_spikes, spike_sum, etc.
% Derived from CC_IV
%
% 1/2001 : Switched spectral analysis to Lomb's PSD method (see Press et al).
%          changed structure to separate traces with and without spikes.
%
% 9/29/03: Cleaned up code and added switches to stop spectral analysis
% if it is not needed. Adjusted array dimensions to match data (previous version was
% transposed, although analysis was still correct).
%
% Paul B. Manis, Ph.D.
% pmanis@med.unc.edu
%


if(nargin > 0)
    source = varargin{1};
    dest = varargin{2};
else
    source = 1;
    dest = 2; % defaults are voltage (src) and current (dst)
end;
%try
sf = getmainselection;
if(sf > 0)
    pflag = getplotflag;
    QueMessage('STA analysis', 1); % clear the que
    for i = 1:length(sf)
        sta_2(sf(i), pflag, source, dest, i);
    end;
end;
%catch
%   watchoff;
%   QueMessage('Error in STA Analysis routine', 1);
%end;


function sta_2(sf, plot_flag, isrc, idst, index)

global DFILE
global ALLCH
global CONTROL

QueMessage('STA - Starting');


dat = [];
time = [];
ivresult = []; % initialize it


% various time definitions:
sta.t_pre = 5; % 35 Time to compute average before a spike; normally 25 msec
sta.t_spec_w = 2; % 5 latest time before spike to measure spectrum (e.g., before foot)
sta.t_post = 20;  % window to display AFTER the spike to show ahp
sta.t_delay = 15; % 20 time after last step before the spike before we consider the data on the rising phase valid
sta.min_isi = 25;
sta.min_isipost = 25;

[DFILE, err] = analysis_setup(DFILE, sf);

if(err ~= 0)
    return;
end;

% flags to control analysis path

sta.do_spike = 1; % enable (=1) detailed analysis of spike shape.
sta.spectgram = 0; % enable (=1) spectrogram analysis (time-dependent spectrum view)
sta.prespikespec = 0; % enable (=1) pre-spike spectrum analysis (windowed spectrum analysis with Lomb's PSD)
sta.nspike = '1:100'; % 10000; % Which spike to go up to in extraction (1 = first, 2 = first and second... etc. -1 means all but first, -2 means all but first and second)
sta.spike_thresh=number_arg(CONTROL(sf).thresh); % spike voltage threshold

if(exist('stapar.mat', 'file') > 0)
    load 'stapar.mat'; % get variables from disk.
end;
if(~ischar(sta.nspike))
    sta.nspike = num2str(sta.nspike);
end;
if(index == 1)
    dlgans = inputdlg({'Pre (ms)', 'Post (ms)', 'Min ISI (ms)', 'Min ISI Post (ms)', 'Thresh', 'Nspike', 'do_spike(0,1)'}, ...
        'STA Parameters', 1, ...
        {num2str(sta.t_pre), num2str(sta.t_post), num2str(sta.min_isi), ...
        num2str(sta.min_isipost), num2str(sta.spike_thresh), sta.nspike, ...
        num2str(sta.do_spike)});
    if(isempty(dlgans))
        return; % probably canceled.
    end;
    sta.t_pre = str2double(dlgans{1});
    sta.t_post = str2double(dlgans{2});
    sta.min_isi = str2double(dlgans{3});
    sta.min_isipost = str2double(dlgans{4});
    sta.nspike = dlgans{6}; % note: nspike can list the spike numbers to be analyzed - [2:5] would do the 2nd to 5th spikes
    sta.spike_thresh = str2double(dlgans{5});
    sta.do_spike = str2double(dlgans{7});
    save 'stapar.mat' sta
end;


QueMessage(sprintf('STA - Spike Threshold = %7.2f mV', sta.spike_thresh));

if(sta.spectgram)
    hmp = findobj('Tag', 'IndividualSpect');
    if(ishandle(hmp))
        hc = get(hmp, 'Children');
        delete(hc); % remove the children plots.
    end;
end;

protocol=deblank(lower(CONTROL(sf).protocol));


% initialize the components of the SPIKE structure
spike.timestamp = datestr(now);
spike.edit = 0; % not editable...
spike.spike_thresh = sta.spike_thresh;
spike.nspike = sta.nspike;
spike.n = 0;
spike.measthresh = NaN; % mean(spikesumm.measthresh); % fraction of height from rest to spike peak where measurement was made
spike.vthresh = NaN; % mean(spikesumm.vthr); % true threshold (2nd derivative method)
spike.vthrmin = NaN; % min(spikesumm.vthr); % min value (first spike...)
spike.deadwin = NaN; % mean(spikesumm.deadwin);
spike.minisi = NaN; % mean(spikesumm.minisi);
spike.minisipost = NaN;
spike.max= NaN; % mean(spikesumm.vmax); spike.max_sd = std(spikesumm.vmax);
spike.apht = NaN; % mean(spikesumm.apht); spike.apht_sd = std(spikesumm.apht);
spike.half_ht = NaN; % mean(spikesumm.half_ht); spike.half_ht_sd = std(spikesumm.half_ht);
spike.min=NaN; % mean(spikesumm.vmin); spike.min_sd = std(spikesumm.vmin);
spike.dvmax = NaN; % mean(spikesumm.dvmax); spike.dvmax_sd = std(spikesumm.dvmax);
spike.dvmin = NaN; % mean(spikesumm.dvmin); spike.dvmin_sd = std(spikesumm.dvmin);
%      dvvmax = spikesumm.dvvmax(find(~isnan(spikesumm.dvvmax)));
spike.dvvmax = NaN; % mean(dvvmax); spike.dvvmax_sd = std(dvvmax);
%      dvvmin = spikesumm.dvvmin(find(~isnan(spikesumm.dvvmin)));
spike.dvvmin = NaN; % mean(dvvmin); spike.dvvmin_sd = std(dvvmin);
%     hwa = spikesumm.hwa(find(~isnan(spikesumm.hwa))); hwb = spikesumm.hwb(find(~isnan(spikesumm.hwb)));
spike.hwa = NaN; % mean(hwa); spike.hwa_sd = std(hwa); % eliminate the ones where coulnd't get a measurement
spike.hwb = NaN; % mean(hwb); spike.hwb_sd = std(hwb);
spike.ahpdepth =  NaN; % mean(spikesumm.ahpdepth); spike.ahpdepth_sd = std(spikesumm.ahpdepth);
spike.fsl = NaN; % first spike latency list...
spike.latency = NaN; % spike latency list
spike.rec = NaN; % spike record list...
spike.tavg = NaN;
spike.vavg = NaN;
spike.vsem = NaN;
spike.spect = NaN;
spike.fspect = NaN;
spike.spect_sem = NaN;
spike.spect_conf = NaN;
spike.pre_record = NaN; % info on time before spike.
spike.pre_slope = NaN;
spike.pre_var = NaN;
spike.pre_specamp = NaN;
spike.pre_freq = NaN;
spike.pre_energy = NaN;

[records, pts]=size(ALLCH{idst});

durs = [number_arg(CONTROL(sf).durho) number_arg(CONTROL(sf).durs1) number_arg(CONTROL(sf).durs2) number_arg(CONTROL(sf).durs3)];

% now compute the minimum spike latency we will consider - this is protocol dependent.

% selections of spikes to analyze
exclb = 0; % 75;
excle = 0; % 140;
exclf = 0; % flag: 0 means no exculsion. % default is not exclusion...

% set latencies for measurements according to the protocol we are using

switch protocol
    case 'forward'
        min_lat = sum(durs(1:3));
        max_dur = sum(durs(1:4));
        % selections of spikes to analyze
        exclb = durs(1); %75;
        excle = sum(durs(1:3)); % 140;
        exclf = 1; % flag: 0 means no exculsion.
    case 'forward2'
        min_lat = 1310; % sum(durs(1:3));
        max_dur = sum(durs(1:4));
        % selections of spikes to analyze
        exclb = 0; %75;
        excle = 1310; % 140;
        exclf = 1; % flag: 0 means no exculsion.

    case {'hyp', 'hypwatch', 'hypsave', 'hypdrug', 'ap-hyp2', 'ap-hyp3'}
        min_lat = sum(durs(1:2));
        max_dur = sum(durs(1:3));

    case {'ap_hypsynch'}
        min_lat = sum(durs(1:2));
        max_dur = sum(durs(1:4));
    case {'drugtestiv'}
        min_lat = 210;
        max_lat = 710;
    otherwise
        min_lat = durs(1); % just the holding time is enough
        max_dur = sum(durs(1:4));
end;

% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
RATES = DFILE.rate .* DFILE.nr_channel / 1000; % array of sampling rates, convert to msec

t0=0;

% get spike times
QueMessage('STA - Finding spikes');
[first_spike, first_isi, nr_spikes, spike_train]=find_spikes2(DFILE, ALLCH{isrc}, t0, tmax, sta.spike_thresh);
QueMessage(sprintf('STA - found %4d spikes', sum(nr_spikes)));
first_isi=max(first_isi,0);

fsl = first_spike(first_spike>0);
fisi = first_isi(first_isi>0);

spike.fsl = fsl;
mean_fsl = mean(fsl); % get mean first spike latency from first spikes in dataset
if(isempty(mean_fsl) || isnan(mean_fsl))
    % try to use spike latencies in data set one above us (assumes that one has spikes)
    mean_fsl = mean(CONTROL(sf-1).spike.fsl);
    if(isempty(mean_fsl) || isnan(mean_fsl)) % still not there?
        mean_fsl = min_lat + sta.t_delay + sta.t_pre; % use a fixed value that is probably the shortest possible...
    end;
end;

pmax=floor((sta.t_pre + sta.t_post)/RATES(1))+1; % this is the max size of the array we need.

tap = {};
vap = {};
iap = {};


% compute the exclusion - don't look at spikes in this window
% note: we do this by modifying the spike list itself... ;)
% only spikes in a given time window.

if(exclf) 
    for i = 1:length(spike_train)
        u = find(spike_train(i).latency < exclb | spike_train(i).latency > excle);
        spike_train(i).latency = [spike_train(i).latency(u)];
    end;
end;
%
% select nth spike in train for analysis 
%  Note limited set of protocols that can use this.
% sta.nspike determines what is done:
% if 0, all spikes are analyzed
% if 1, it analyzes only the first spike
% if 2, it analyzes only the second spike
% if [1:12], it analyzes the first throug the 12th spikes...
% if [2:1000] it analyzes all but the first spike, 
% if [1:1000] it analyzes all the spikes
% etc. 
% but can be second 
if(strmatch(CONTROL(sf).protocol,{'hyp2', 'forward2', 'hyp_iv'}, 'exact')) % then only compute averages for the Nth spike
    stim_list = record_parse(CONTROL(sf).sequence);
    nstm = length(stim_list);
    rep_length = floor(records/nstm);
    QueMessage(sprintf('STA - seq: %s  nstm: %d  reps: %d', CONTROL(sf).sequence, nstm, rep_length));
%    fprintf(1, 'Spike trains: %d\n', length(spike_train));
    for i = 1:length(spike_train)
        thisNSpikes = length(spike_train(i).latency);
%        fprintf(1, '   Train %d has %d spikes (to start), sta.nsp = %s\n', i, thisNSpikes, sta.nspike );
        r = spike_train(i).source;
        istim = 1 + floor((r-1)/rep_length);
        if(istim <= length(stim_list))
            maxspike = min(max(str2num(sta.nspike)), thisNSpikes); % spike selection code - top end
            minspike = min(str2num(sta.nspike)); % earliest by count
            u = find(spike_train(i).latency(minspike:end) > min_lat)+minspike-1;  %#ok<MXFND>
            splist = min(u):min(min(max(u), maxspike), length(spike_train(i).latency)); % make sure it is in the list
%            fprintf(1, '    max spike: %d    selected spikes %s (#=%d)\n', maxspike, num2str(splist), length(splist));
            spike_train(i).latency = [spike_train(i).latency(splist)];
        else
        spike_train(i).latency = [];
        end;
    end;

end;

% copy results to ISILAT
v = size(spike_train);
ISILAT(1:v(2))=deal(struct('nspikes', 0, 'isi', 0, 'lat', 0, 'source' , 0));
% initialize the ISILAT structure
for i=1:v(2)
    ISILAT(i).nspikes = length(spike_train(i).latency);
    ISILAT(i).isi = diff(spike_train(i).latency);
    ISILAT(i).lat = spike_train(i).latency;
    ISILAT(i).source = spike_train(i).source; % save the source trace number too
end


% setup for spectrum analysis

nfft =2048;
Pvgrand = [];
Psgrand = [];
flomb = [1:300];
Plgrand = zeros(records, length(flomb));
Pv=[];
Pvl=[];
ngrand = 0;
conf_t = 0; % average confidence limit
prob = 0.05; % use 5% level to set


if(sta.prespikespec)
    QueMessage('Building multiplot');
    % just for testing to see what kind of data is being analyzed..
    fsize = 7;
    ng = ceil(sqrt(size(ALLCH{isrc}, 1)));
    hmp = multiplot('IndividualSpect', ng, ng, [0.6, 0.4]);
end;

% **************************
QueMessage('Starting analysis loop');
k = 0;
v_npoints = size(ALLCH{isrc}, 2);

% Initialize the output array information

fsl = NaN*zeros(records, 1);
nospk = zeros(records, 1);
pre_vm = fsl;
pre_var = fsl;
pre_maxmv = fsl;
pre_slope = fsl;
pre_freq = fsl;
pre_specamp = fsl;
pre_energy = fsl;
pre_record = [1:records];
acceptn = 0;
for rec = 1:records
    spl = [];
    if(exist('spike_train', 'var') > 0 && ~isempty(spike_train) && any(strcmp('source', fieldnames(spike_train))))
        spl = find([spike_train.source] == rec); % match spikes in the list
    end;
    if(~isempty(spl)) % if there are spikes
        for i=[spl]  % for every spike in the list
            for j = 1:length(spike_train(i).latency)  % check all spikes in the response
                doit = 0; % flag to determine whether we analyze this spike or not
                if(spike_train(i).latency(j) > (min_lat + sta.t_pre + sta.t_delay))
                    if(j == 1)
                        doit = 1; % always do the first spike if it meets the criteria
                    elseif((spike_train(i).latency(j) - spike_train(i).latency(j-1)) > sta.t_pre)
                        doit = 1; % disallow preceding spikes within the averaging window
                    end;
                end;
                if(doit) % analysis on the spike
                    if(isnan(fsl(rec))) % have we identified a first spike yet?
                        fsl(rec) = spike_train(i).latency(j) - min_lat; % reference to start of step
                    end;
                    nospk(rec) = NaN;
                    pb0 = (spike_train(i).latency(j) - sta.t_pre)/RATES(rec); % compute start of spike
                    pb = floor(pb0)+1;
                    dtfix = RATES(rec)*(pb-pb0); % calculate offset for this trace - exact alignment on threshold crossing
                    pe = floor((spike_train(i).latency(j) + sta.t_post)/RATES(rec))+1; % compute end of spike window
                    if(((pe-pb+1) - pmax) > 0 && pe < v_npoints && pb > 0) % only work with spikes where complete data is available.
                        k = k + 1;
                        vap{k} = ALLCH{isrc}(rec, pb:pe); % cut out the spike from the window
                        npts = length(vap{k});
                        tap{k} = (([0:npts-1]*RATES(rec)) - sta.t_pre + dtfix); % offset time by zero crossing amount (interpolated)
                        iap{k} = ALLCH{idst}(rec, pb:pe);
                        pw = floor((spike_train(i).latency(j) - sta.t_spec_w)/RATES(rec))-pb; % compute pre-spike window

                        V = vap{k}(1:pw);
                        T = tap{k}(1:pw);
                        % estimate rising slope and noise (variance) during the slope
                        p1 = polyfit(T, V, 1); % fit a line also
                        pre_slope(rec) = p1(1); % store slope value
                        pre_vm(rec) = mean(V);
                        xl = V - (p1(1)*T(1) + p1(2)); % linear subtraction

                        pre_var(rec) = sqrt(var(xl)); % compute variance around the mean rise
                        pre_maxmv(rec) = max(abs(xl)); % max deviation from the line
                        pre_record(rec) = rec;
                        ngrand = ngrand + 1;

                        if(sta.prespikespec) % computations for pre-spike spectrum
                            % compute power spectrum of this phase:
                            % first subtract baseline as a second - order polynomial
                            p = polyfit(T, V, 2);
                            Vp = V - (p(1)*T.^2 +p(2)*T+ p(3));
                            % window the data.
                            w = hanning(length(Vp), 'symmetric');
                            win_mag = sum(w)/length(w);
                            Vp = Vp .* w;
                            n_pts = floor(length(Vp)/2);
                            % Lomb's PSD method
                            [fl, Pvl, conf] = lombpsd2(T/1000, Vp, flomb, prob);
                            Pvl = sqrt((Pvl * 2 * var(Vp))/n_pts)/win_mag; % remove normalization and correct for window
                            conf = conf / (RATES(rec)*pw);
                            Plgrand(rec,:) =  Pvl;
                            conf_t = conf_t + conf;
                            [pre_specamp(rec), ipsdm] = max(Pvl); % store peak amplitude and frequency
                            pre_freq(rec) = fl(ipsdm);
                            pre_energy(rec) = sum (Pvl); % total energy is sum of spectral content over frequency...
                        end;


                        if(sta.spectgram) % computations for spectrogram prior to spike
                            axes(hmp(rec, 1));
                            plot(fl, Pvl, '-k');
                            if(mod(rec, ng) ~= 1)
                                xlabel('');
                                set(gca, 'XTickLabel', '');
                            end;
                            axes(hmp(rec, 2));
                            plot(T/1000, Vp, '-k'); % plot what is fit to psd
                            hold on;
                            plot(T/1000, xl, '-r');
                            li = p1(1)*T + p1(2);
                            set(gca, 'Fontsize', 6, 'XTickLabel', '');
                            if(mod(rec, ng) ~= 1)
                                ylabel('');
                                set(gca, 'YTickLabel', '');
                            end;

                            % time-spectrogram
                            [Ps, fs] = specgram(ALLCH{isrc}(k,:), 256, fsamp, 4, 3);
                            if(isempty(Psgrand))
                                Psgrand = zeros(size(Ps,1),size(Ps,2));
                            end;
                            Psgrand = Ps + Psgrand;
                        end;
                    end;
                end;
            end;
        end;
    else % no spikes in the record - calculate slope, variance and spectrum anyway...
        %	insert "problematic code" here.....
    end;
end;
QueMessage('Completed analysis loop');

spmax = 0;
vmax = 0;
if(sta.spectgram)
    for i = 1:length(pre_record)
        axes(hmp(i, 1));
        u = get(gca, 'YLim');
        if(u(2) > spmax)
            spmax = u(2);
        end;
        axes(hmp(i, 2));
        u = get(gca, 'Ylim');
        if(max(abs(u)) > vmax)
            vmax = max(abs(u));
        end;
    end;
    for i = 1:length(pre_record)
        axes(hmp(i, 1));
        set(gca, 'YLim', [0, spmax]);
        axes(hmp(i, 2));
        set(gca, 'YLim', [-vmax vmax]);
        text(0, 0.7*vmax, sprintf(' %d', pre_record(i)), 'FontSize', 7, 'Color', 'red');
    end;
end;
Plave = [];
Plsem = [];
conf_ave = [];

if(ngrand > 0)
    %   Pvave = mean(Pvgrand'); % /ngrand;
    %   Pvsem = std(Pvgrand')/sqrt(size(Pvgrand, 1));
    %   [Pfpat, Pvpat] = make_patch(Pvf(:,2), Pvave, Pvsem);
    if(sta.prespikespec)
        Plave = mean(Plgrand);
        Plsem = std(Plgrand)/sqrt(size(Plgrand, 1));
        [Plfpat, Plvpat] = make_patch(flomb, Plave, Plsem);
        conf_ave = conf_t / ngrand; % average confidence level
    end;
    if(sta.spectgram)
        Psgrand = Psgrand/ngrand; % get average
    end;

    % flush out the ones where we didn't have any spikes with an appropriate window
    % interpolate all spikes onto a standard first time base
    QueMessage ('Interpolating');
    tap_ref = [0:npts-1]*RATES(1) - sta.t_pre;
    n=size(tap, 2);
    for i = 1:n
        vapi(i,:) = interp1([tap{i}], [vap{i}], tap_ref);
        iapi(i,:) = interp1([tap{i}], [iap{i}], tap_ref);
    end;
    %   vap=[];
    % iap=[];
    nr=length(tap_ref);

    vapi = vapi(:,2:end-1);
    iapi = iapi(:,2:end-1);
    tapn=tap_ref(2:end-1);
    QueMessage('Averaging');
    nspikes=size(vapi, 1);
    vavg = mean(vapi);
    vstd = std(vapi)/sqrt(nspikes); % converted to SEM 10/2/2000 P. Manis
    iavg = mean(iapi);
    istd = std(iapi)/sqrt(nspikes);
    tavg = tapn;

    % calculate the mean current variation based on Bryant and Segundo
    % (1976). This is 2 SE's of the mean of the current variance for all
    % traces. Because the current is ramped up and doesn't reach steady
    % state until about 100 msec, we use the window 100-900 msec for this
    % calculation.

    t0=floor(100/RATES(1));
    t1=floor(400/RATES(1));
    if(t1 > size(ALLCH{1}, 2))
        t1 = size(ALLCH{1}, 2);
    end;
    inj_mean = mean(mean(ALLCH{idst}(:,t0:t1))); % that is the true mean current (all traces)
    c=reshape(ALLCH{idst}(:,t0:t1), [1,numel(ALLCH{idst}(:,t0:t1))]);
    inj_conflim = 2*std(c)/sqrt(size(vapi,1));
    clear c  % remove extra copy of current.
    %    fprintf(1, 'Imean: %9.3f   Iconflim: %9.4f\n', inj_mean, inj_conflim);
    %    inj_mean
    %   inj_conflim

    [tpat, vpat] = make_patch(tavg, vavg, vstd);
    [tpat, ipat] = make_patch(tavg, iavg, istd);



    % analyze the spikes - hw, dv/dt, amplitude from rest, etc.

    spike_count = size(ISILAT);
    if(spike_count(2) > 0 && sta.do_spike == 1)
        QueMessage('STA - Summarizing Spikes');
        spikesumm = spike_sum(DFILE, ALLCH{isrc}, ISILAT, 0, CONTROL(sf).Rmp, CONTROL(sf).filename, CONTROL(sf).recbeg, CONTROL(sf).deadwin, ...
            [sta.min_isi sta.min_isipost], 3, 10, plot_flag);

        % summary spike information:
        if(~isempty(spikesumm))
            n = size(spikesumm.vmax);
            spike.timestamp = datestr(now);
            spike.edit = 0; % not editable...
            spike.n = n(2);
            spike.measthresh = mean(spikesumm.measthresh); % fraction of height from rest to spike peak where measurement was made
            athr = find(~isnan(spikesumm.vthr));
            if(isempty(athr))
                spike.vthresh=NaN;
                spike.vthrmin=NaN;
            else
                spike.vthresh = mean(spikesumm.vthr(athr)); % true threshold (2nd derivative method)
                spike.vthrmin = min(spikesumm.vthr(athr)); % min value (usually the first spike...)
            end;
            spike.n = length(spikesumm.vmax);
            spike.deadwin = mean(spikesumm.deadwin);
            spike.minisi = mean(spikesumm.minisi);
            spike.minisipost = mean(spikesumm.minisipost);
            spike.max=mean(spikesumm.vmax); spike.max_sd = std(spikesumm.vmax);
            spike.apht = mean(spikesumm.apht); spike.apht_sd = std(spikesumm.apht);
            spike.half_ht = mean(spikesumm.half_ht); spike.half_ht_sd = std(spikesumm.half_ht);
            spike.min=mean(spikesumm.vmin); spike.min_sd = std(spikesumm.vmin);
            spike.dvmax = mean(spikesumm.dvmax); spike.dvmax_sd = std(spikesumm.dvmax);
            spike.dvmin = mean(spikesumm.dvmin); spike.dvmin_sd = std(spikesumm.dvmin);
            dvvmax = spikesumm.dvvmax(find(~isnan(spikesumm.dvvmax)));
            spike.dvvmax = mean(dvvmax); spike.dvvmax_sd = std(dvvmax);
            dvvmin = spikesumm.dvvmin(find(~isnan(spikesumm.dvvmin)));
            spike.dvvmin = mean(dvvmin); spike.dvvmin_sd = std(dvvmin);
            hwa = spikesumm.hwa(find(~isnan(spikesumm.hwa))); hwb = spikesumm.hwb(find(~isnan(spikesumm.hwb)));
            spike.hwa = mean(hwa); spike.hwa_sd = std(hwa); % eliminate the ones where coulnd't get a measurement
            spike.hwb = mean(hwb); spike.hwb_sd = std(hwb);
            spike.ahpdepth =  mean(spikesumm.ahpdepth); spike.ahpdepth_sd = std(spikesumm.ahpdepth);


        end;
    end;
    QueMessage('STA - Storing Spike Info');
    spike.tavg = tavg;
    spike.vavg = vavg;
    spike.vsem = vstd; % actually, SEM.
    spike.iavg = iavg;
    spike.isem = istd;
    spike.injconf = inj_conflim;
    spike.injmean = inj_mean;
    spike.fsl = fsl;
    spike.latency = {spike_train.latency};
    spike.source = {spike_train.source};
    spike.spect = Plave; % prespike spectrum (using Lomb's method for 0-300 Hz.)s
    spike.spect_sem = Plsem; % sem of pre spike spectrum (across records)
    spike.spect_conf = conf_ave;
    spike.fspect = flomb; % and frequency plot for that
    spike.pre_record = pre_record;
    spike.pre_vm = pre_vm;
    spike.pre_slope = pre_slope;
    spike.pre_var = pre_var;
    spike.pre_freq = pre_freq;
    spike.pre_specamp = pre_specamp;
    spike.pre_energy = pre_energy;
spike.spike_thresh = sta.spike_thresh;
spike.nspike = sta.nspike;
spike.ngrand = ngrand;

    if(sta.prespikespec)
        hs = findobj('tag', 'V-S');
        if(isempty(hs) | ~ishandle(hs))
            hs = figure('Tag', 'V-S', 'name', 'vm vs spec amp', 'NumberTitle', 'off');
        else
            figure(hs);
        end;
        subplot(2,1,2);
        plot(spike.pre_vm, spike.pre_specamp,'rx');
        subplot(2,1,1);
        plot(spike.pre_vm, spike.pre_slope, 'bo');
    end;

    CONTROL(sf).spike=spike; % make it a substructure
else
    QueMessage('STA: ? no spikes meet criteria');
    return;
end;
pc=0.5; % patch color (grey)

% spike train stuff...
if(plot_flag >= 0)
    if(sta.prespikespec)
        h = findobj('Tag', 'STA_spec');
        if(isempty(h))
            h = figure('Tag', 'STA_spec', 'name', 'pre-spike spectrum', 'NumberTitle', 'off');
        end;
        figure(h);
        clf;
        fsize = 7;
        subplot('Position', [0.07, 0.55, 0.43, 0.4]);
        patch(Plfpat, Plvpat, [pc pc pc], 'EdgeColor', 'none');
        hold on;
        plot(flomb, Plave, '-b');
        hold on;
        plot([fl(1) fl(end)], [conf_ave conf_ave], 'k--');
        %   set(gca, 'XScale', 'log');
        [m,k]=max(Plave);
        ylabel('mV/Hz^{1/2}', 'fontsize', fsize);
        xlabel('F (Hz)', 'fontsize', fsize);
        text(fl(k), Plgrand(k), sprintf('\\uparrow\nFo = %8.1f Hz\nPer = %8.1f ms\nMax = %8.2f mV\^2/Hz', ...
            fl(k), 1000/fl(k), sqrt(Plave(k))),...
            'FontSize', 6, 'VerticalAlignment', 'Top', 'HorizontalAlignment', 'Center');
        % fill in other squares... with record vs. peak f, energy measure, slope and variance.

        subplot('Position', [0.57, 0.77, 0.43, 0.2]);
        plot(pre_record, pre_freq, 'ko', 'MarkerFaceColor', 'black', 'MarkerSize', 2.5);
        ylabel('Peak F (Hz)', 'fontsize', fsize);
        subplot('Position', [0.57, 0.57, 0.43, 0.2]);
        plot(pre_record, pre_specamp, 'ro', 'MarkerFaceColor', 'red', 'MarkerSize', 2);
        ylabel('Amplitude', 'fontsize', fsize);
        xlabel('Record', 'fontsize', fsize);

        subplot('Position', [0.07, 0.32, 0.43, 0.15]);
        plot(pre_record, fsl, 'ko', 'MarkerFaceColor', 'black', 'MarkerSize', 2.5);
        hold on;
        plot(pre_record, nospk, 'ro', 'MarkerFaceColor', 'red', 'MarkerSize', 2);
        ylabel('fsl (msec)', 'fontsize', fsize);

        subplot('Position', [0.07, 0.07, 0.43, 0.25]);
        plot(pre_record, pre_slope, 'ko', 'MarkerFaceColor', 'black', 'MarkerSize', 2.5);
        u=get(gca, 'Xlim');
        set(gca, 'Xlim', [0 u(2)]); % set min to 0
        xlabel ('Record', 'fontsize', fsize);
        ylabel ('Slope (mV/ms)', 'fontsize', fsize);

        subplot('Position', [0.57, 0.27, 0.43, 0.2]);
        plot(pre_record, pre_var, 'ko', 'MarkerFaceColor', 'black', 'MarkerSize', 2);
        hold on;
        plot(pre_record, pre_maxmv, 'b^', 'MarkerFaceColor', 'blue', 'MarkerSize', 3);
        ylabel('Variance (mV), max deviation', 'fontsize', fsize);

        subplot('Position', [0.57, 0.07, 0.43, 0.2]);
        plot(pre_record, pre_energy, 'ro', 'MarkerFaceColor', 'red', 'MarkerSize', 2);
        ylabel('Energy', 'fontsize', fsize);
        xlabel('Record', 'fontsize', fsize);
    end;

    if(sta.spectgram)
        h = findobj('Tag', 'STA_specgm');
        if(isempty(h))
            h = figure('Tag', 'STA_specgm', 'name', 'Spectrogram', 'NumberTitle', 'off');
        end;
        figure(h);
        clf;
        imagesc(20*log10(abs(Psgrand))); axes xy; colormap(jet);
    end;

    % now the main figure
    h = findobj('Tag', 'STA_spiketrain'); % check for pre-existing window
    if(isempty(h)) % if none, make one
        h = figure('Tag', 'STA_spiketrain', 'Name', 'spike train Analysis', 'NumberTitle', 'off');
    end
    figure(h); % otherwise, select it
    clf; % always clear the window...
    fsize = 7;
    msize = 3;
    subplot('Position', [0.1 0.55 0.35 0.35]);

    stl = 0;
    if(exist('spike_train') & ~isempty(spike_train) & any(strcmp('latency', fieldnames(spike_train))))
        stl = 1;
        sll = []; trial=[];

        for i = 1:length(spike_train)
            sll = [sll spike_train(i).latency'];
            u = length(spike_train(i).latency); % how many spikes?
            trial = [trial [zeros(u,1)+i]'];
        end;
        plot(sll, trial, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 2.7);
    end;
    u = get(gca, 'Xlim');
    u(1)=0;
    set(gca, 'Xlim', u);
    subplot('Position', [0.1 0.1 0.35 0.35]);
    if(stl)
        hist(sll, 50);
        set(gca, 'Xlim', u);
    end;
    subplot('Position', [0.55 0.55 0.35 0.35]);
    %         plot(tap, vap);
    if(strmatch(CONTROL(sf).protocol,{'hyp2', 'hyp', 'hypdrug'}))
        %      plot(tavg, vapi, 'linewidth', 0.1);
    else
        %      plot(tavg, iavg);
    end;
    pc=0.85;
    patch(tpat, ipat-inj_mean, [pc pc pc], 'EdgeColor', 'none');
    hold on;
    plot([tavg(1) tavg(end)], [inj_conflim inj_conflim], 'r');
    plot([tavg(1) tavg(end)], [-inj_conflim -inj_conflim], 'r');
    plot(tavg, iavg-inj_mean, 'b');
    grid;

    subplot('Position', [0.55 0.1 0.35 0.35]);
    pc=0.85;
    patch(tpat, vpat, [pc pc pc], 'EdgeColor', 'none');
    hold on;
    plot(tavg, vavg);
    grid;

    subplot('Position',[0.1,0.92,0.8,0.08])
    axis([0,1,0,1])
    axis('off')
    text(0,0.80,sprintf('%-12s R[%d:%d]     %-8s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Fontsize', 8);
    text(0,0.5,sprintf('Gain:%4.1f  LPF:%4.1f kHz   Pre: %6.1f ms  Post: %6.1f ms  Min ISI: %6.1f ms', DFILE.igain, DFILE.low_pass(1), ...
        sta.t_pre, sta.t_post, sta.min_isi), 'FontSize', 8);
    text(0, 0.2, sprintf('# spikes in avg: %d   Mean current: %8.3f pA  conf: %9.3f pA', nspikes, inj_mean, inj_conflim), 'Fontsize', 8);
    text(0.7,0.8,sprintf('%s (sta code: 9/30/03)',date),'FontSize',8);
    orient landscape

    % control printing and closing of window for automatic runs
    % plot_flag = 1 creates plot and leaves it up
    % plot_flag = 2 creates plot but closes it when done
    % plot_flag = 2 creates plot and prints it and then closes it
    if (plot_flag == 2)
        print -dljet3;
    end
    if plot_flag >= 2
        close
    end

end;
return;

function problematic()
if(rogerdidit)
    % this is mostly a repeat of the code above, but we don't have to worry about spikes.
    pb = floor((mean_fsl - sta.t_pre)/RATES(rec)); % absolute time - start this display at min_measure msec.
    pe = floor((mean_fsl + sta.t_post)/RATES(rec)); % match to data above...
    pw = floor((mean_fsl - sta.t_spec_w)/RATES(rec))-pb; % compute pre-spike window
    k = k + 1;
    vap(:,k) = VOLTAGE(rec, pb:pe)'; % cut out the spike from the window
    npts = length(vap);
    tap(:,k) = (([0:npts-1]*RATES(rec)) - sta.t_pre)';
    iap(:,k) = CURRENT(rec, pb:pe)';
    V = vap(1:pw, k);
    T = tap(1:pw, k);
    p1 = polyfit(T, V, 1); % fit a line also
    pre_slope(rec) = p1(1); % store slope value
    pre_vm(rec) = mean(V);
    xl = V - (p1(1)*T + p1(2));
    pre_var(rec) = sqrt(var(xl)); % compute variance around the mean rise
    pre_maxmv(rec) = max(abs(xl)); % max deviation from the line
    pre_record(rec) = rec;
    ngrand = ngrand + 1;

    if(sta.prespikespec)
        % compute power spectrum of this phase:
        % first subtract baseline as a second - order polynomial
        p = polyfit(T, V, 2);
        Vp = V - (p(1)*T.^2 +p(2)*T+ p(3));
        % window the data.
        w = hanning(length(Vp), 'symmetric');
        n_pts = floor(length(Vp) / 2)
        win_mag = sum(w)/length(w);
        Vp = Vp .* w;
        % Lomb's PSD method
        [fl, Pvl, conf] = lombpsd2(T/1000, Vp, flomb, prob); % we measure it, and add to "pre-spike" data,
        Pvl = sqrt(Pvl * (2 * var(Vp)) / n_pts) / win_mag; % remove normalization provided by lombpsd
        conf = conf / (RATES(rec)*length(T));
        Plgrand(rec,:) = Pvl;

        [pre_specamp(rec), ipsdm] = max(Pvl); % store peak amplitude and frequency
        pre_freq(rec) = fl(ipsdm);
        pre_energy(rec) = sum (Pvl); % total energy is sum of spectral content over frequency...
    end;

    if(sta.spectgram)
        axes(hmp(rec, 1));
        plot(fl, sqrt(Pvl), '-k'); %  fl, conf, '--r');
        if(mod(rec, ng) ~= 1)
            xlabel('');
            set(gca, 'XTickLabel', '');
        end;
        if(rec < (ng*(ng-1)))
            ylabel('');
            set(gca, 'YTickLabel', '');
        end;
        axes(hmp(rec, 2));
        plot(T/1000, Vp, '-k');
        hold on;
        li = p1(1)*T + p1(2);
        set(gca, 'Fontsize', 6);
        if(mod(rec, ng) ~= 1)
            xlabel('');
            set(gca, 'XTickLabel', '');
        end;
    end;
end;

